package fem2;

import math2.FunctionBasisOnRnToR;
import math2.FunctionRnToR;
import math2.LinAlgUtilities;
import math2.PolynomialBasisOnRnToR;
import fem2.enu.MeshPartType;

/**
 * Represent a distributed load over a mesh part (edge, face or block)
 * 
 * @author hbui
 * 
 */
public class DistributedLoad extends Load {

	private FunctionRnToR[] lcs;

	/**
	 * 
	 * @param mp
	 * @param lcs
	 */
	public DistributedLoad(MeshPart mp, FunctionRnToR... lcs) {
		super(mp);
		assert lcs.length == mp.countFields() : "incompatible load components and active fields on mesh part";
		this.lcs = lcs;
	}

	@Override
	public double[] computeRe(double t) {
		MeshPart mp = getMeshPart();

		FunctionBasisOnRnToR basis = mp.getBasis();
		int n = basis.getDimension();
		int d = mp.countFields();
		double[] re = new double[n * d];

		/*
		 * select the maximum order for integration
		 */
		int p;
		if (basis instanceof PolynomialBasisOnRnToR) {
			p = ((PolynomialBasisOnRnToR) basis).getP();
		} else {
			p = 5; // set a large p to account for accuracy
		}
		int pn = 0;
		for (int j = 0; j < d; j++) {
			if (this.lcs[j] != null) {
				if (lcs[j].getP() > pn) {
					pn = lcs[j].getP();
				}
			}
		}
		p += pn;

		/*
		 * integration using Gauss quadrature
		 */
		IntegrationManager integrator = IntegrationManager.getInstance();
		MeshPartType type = mp.getType();
		for (int i = 0; i < n; i++) {
			FunctionRnToR Ni = basis.getBasisFunction(i);
			for (int j = 0; j < d; j++) {
				if (this.lcs[j] != null) {
					int ng = integrator.getNumIntegrationPoints(type, p);
					for (int k = 0; k < ng; k++) {
						double[] xi = integrator.getIntegrationPoint(type, p, k);
						double w = integrator.getIntegrationWeight(type, p, k);
						double[][] jac = mp.jacXAt(xi);
						double detJac = Math.sqrt(LinAlgUtilities.determinantJJt(jac));
						re[i * d + j] += lcs[j].valueAt(xi) * Ni.valueAt(xi) * (w * detJac)
								* lambda.valueAt(t);
					}
				} else {
					re[i * d + j] = 0;
				}
			}
		}
		return re;
	}
}
